%==============================================================================
% This code computes equilibrium with CBDC under different CBDC rate
%  COPY RIGHT: Yu Zhu
%==============================================================================


clear;close all
options = optimset('maxiter',5e5,'tolX',1e-10,'tolfun',1e-10,'maxfunevals',5e5);

%==========================================
%Define parameters and functions
%==========================================
name = 'LN87_08_final' 
load(['pars_post_' name '.mat'])
location='figures/'

alpha1=pars(4);alpha2=pars(5);alpha3=pars(6);beta=pars(7);epsilon=0.001;sigma=pars(1);
cost=pars(8);
B=pars(2);
A=pars(10);  %productivity
eta=pars(9);                                       %bank's cost to manage deposits
theta=pars(11);
rq_ratio=pars(12); %reserve ratio
i_reserve = pars(13);
year_option = pars(14); % 1: using banking data from 2007 and 2008. 0: using banking data from 2014-2019
%% Define the name of the saved files: last two number indicates reserve requirement ratio and the interest on reserves
if year_option==1
    pi_m = 0.03346;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100)),'_', num2str(beta*100), '_07_08'];
else
    pi_m = 0.01515;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100)),'_', num2str(beta*100)];
end




                                                    %default rate
Nb=pars(3); %number of banks
nrhogrid =10000;
i_m = 1/beta*(1+pi_m)-1;i_h = 1/beta*(1+pi_h)-1;
pi_eff = (1+pi_m)/(1+i_reserve)-1;

u = @(x) ((x+epsilon).^(1-sigma)-epsilon.^(1-sigma))./(1-sigma); %DM utility
du = @(x)(x+epsilon).^(-sigma); %derivative of DM utility
d2u=@(x)(x+epsilon).^(-sigma-1)*(-sigma);
c=@(x)x;                                                         % DM cost
dc=@(x)1;%   derivative of DM cost
d2c=@(x)0;

y_star = fminsearch(@(x)(du(x)-dc(x)).^2,1);                     % y_star
g_func=@(x) (1-theta).*u(x)+theta.*c(x);                         %trading proctocol 
dg_func=@(x) (1-theta).*du(x)+theta.*dc(x);
d2g_func=@(x)(1-theta).*d2u(x)+theta.*d2c(x);
y_grid = 0:0.001:y_star;
p_grid = g_func(y_grid);
p_star = p_grid(end);
y_func = @(p)interp1(p_grid,y_grid,min(p,max(p_grid)));          %trading quantity 
p_func=@(x)min(x,max(p_grid));                                   %trading payment
lambda_func = @(z) max(du(y_func(z))./dg_func(y_func(z))-1,0);          % liquidity premium

f=@(k)A*k.^(eta);%production
df=@(k)A*eta*k.^(eta-1); %marginal production
df_inv = @(rho)(A*eta./(1+rho)).^(1/(1-eta));


dlambda_func=@(z)((d2u(y_func(z)).*dg_func(y_func(z))-du(y_func(z)).*d2g_func(y_func(z)))./dg_func(y_func(z)).^3).*(z<=p_grid(end));
dinv=@(z)beta*(alpha2*lambda_func(z)+1);
ddinv=@(z)beta*alpha2*dlambda_func(z);
z_grid = linspace(0,p_grid(end),100000);
lambda_grid = lambda_func(z_grid);
gamma_grid=linspace(0,20,nrhogrid);
x0=[0,0,0];
x0_R=x0;
opt=1;
icbdc_grid=linspace(-0.01,(1+pi_m)/beta-1,1000);
picbdc=pi_m;
D_grid = linspace(0,max(p_grid),5e3);
pass_variables
ygrid = D_inv_demand(i_m,D_grid,func_var);
rho_grid = linspace(-.05,rho_bar,1000);
[LS,LD,rho_grid] = LSLD(rho_grid,func_var,pi_m,Nb,ygrid);
LD_cournot = A*eta*(((eta-1)/Nb+1)./(1+rho_grid)).^(1/(1-eta));
eq_plot=figure
hold on
ID = find(LS<=0,1,'last');
h1=plot(1+rho_grid(1:ID+1),[LS(1:ID),0],'black','linewidth',2)
plot(1+rho_grid(ID+1:end),LS(ID+1:end),'black','linewidth',2)
plot([1+rho_grid(ID+1),1+rho_grid(ID+1)],[0,LS(ID+1)],'black','linewidth',2)
plot([1+rho_grid(end),1+rho_grid(end)],[LS(end),LS(end)*2],'black','linewidth',2)
h2=plot(1+rho_grid, LD,'blue','linewidth',2)
%plot(1+rho_grid, LD_cournot,'linewidth',2)

ylim([0,1.7])
xlabel('$R_{\ell}$','interpreter','latex')
ylabel('$L$','interpreter','latex')
legend([h1,h2],'Loan Supply','Loan Demand','location','northwest')
figure_halfpage
saveas(eq_plot,[location,'eq_plot'],'epsc')                                 % generate Figure 5b in the paper
saveas(eq_plot,[location, 'eq_plot'],'png')

for i=1:length(icbdc_grid)
    i
    
    if i==1
    y=SS_eq_noCBDC(func_var,pi_eff,Nb,ygrid);
    end
    eq_picbdc=(1+picbdc)/(1+icbdc_grid(i))-1;
    y_reserve = SS_eq_CBDCR(func_var,pi_m, i_reserve, icbdc_grid(i),Nb,ygrid);
    
    i_m=(1+pi_m)/beta-1;
    psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
    drate(i)=1/psi-1;
    loan= fzero(@(x)df(x)-1-y(end),[1e-10,100]);
    
    eq_final_nocbdc(i,:)=[y,drate(i),loan];
    x0=y;
    x0_R=y_reserve;
    psi_reserve =  beta*(alpha2*lambda_func(y_reserve(2))+alpha3*lambda_func(y_reserve(1)+y_reserve(2)))+beta;
    drate_reserve(i)=1/psi_reserve-1;
    loan_reserve=fzero(@(x)df(x)-1-y_reserve(end),[1e-10,100]);
    
    %===============================================
    %Compute equilibrium where CBDC is not a reserve
    %===============================================
    icbdc_real  = (1+picbdc)/(icbdc_grid(i)+1)/beta-1;
    if icbdc_real>psi/beta-1
        eq_cbdc(i,:)=[y,drate(i),loan];
        drate_cbdc(i)=drate(i);
        eliquidity(i,1) = y(2);
    else
        tic
        z_eq= SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,i_m,icbdc_real);
        
        lsupply_max=(1-rq_ratio)*(1+icbdc_real)*beta*z_eq(1,2);
        rho_hat=max(1/(1+pi_eff)-1, ((icbdc_grid(i)+1)/(1+picbdc)-rq_ratio/(1+pi_eff)+cost)/(1-rq_ratio)-1);
        ld_low= fzero(@(x)df(x)-1-rho_hat,[1e-10,100]);
        drate_cbdc(i)=(icbdc_grid(i)+1)/(1+picbdc)-1;
        toc
        if ld_low<lsupply_max
            deposit =  ld_low/(icbdc_real+1)/beta/(1-rq_ratio);
            eq_cbdc(i,:)=[z_eq(1),deposit,rho_hat,drate_cbdc(i),ld_low];
        else
            rhocbdc = df(lsupply_max)-1;
             
            eq_cbdc(i,:)=[z_eq,rhocbdc,drate_cbdc(i),lsupply_max];
            
        end
        eliquidity(i,1) = z_eq(2);

    end
    
    %===============================================
    %Compute equilibrium where CBDC is a reserve
    %===============================================
    if icbdc_real>psi_reserve/beta-1
        eq_cbdc_R(i,:)=[y_reserve,drate_reserve(i),loan_reserve];
        drate_cbdc_R(i)=drate_reserve(i);
        eliquidityR(i,1) = y_reserve(2);
    else
        tic
        z_eq= SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,i_m,icbdc_real);
        
        lsupply_maxR=(1-rq_ratio)*(1+icbdc_real)*beta*z_eq(1,2);
        rho_hatR=max(1/(1+pi_eff)-1, ((icbdc_grid(i)+1)/(1+picbdc)-rq_ratio*max(1/(1+pi_eff),1/(1+eq_picbdc))+cost)/(1-rq_ratio)-1);
        ld_lowR= fzero(@(x)df(x)-1-rho_hatR,[1e-10,100]);
        drate_cbdc_R(i)=(icbdc_grid(i)+1)/(1+picbdc)-1;
        toc
        if ld_lowR<lsupply_maxR
            depositR =  ld_lowR/(icbdc_real+1)/beta/(1-rq_ratio);
            eq_cbdc_R(i,:)=[z_eq(1),depositR,rho_hatR,drate_cbdc_R(i),ld_lowR];

            %break
        else
            rhocbdcR = df(lsupply_maxR)-1;
             
            eq_cbdc_R(i,:)=[z_eq,rhocbdcR,drate_cbdc_R(i),lsupply_maxR];
            
        end
        eliquidityR(i,1) = z_eq(2);

    end
    
end
markup=(alpha1*min(eq_cbdc(1,1),max(p_grid))/y_func(eq_cbdc(1,1)) +alpha2*min(eq_cbdc(1,2),max(p_grid))/y_func(eq_cbdc(1,2))...
    +alpha3*min(eq_cbdc(1,2)+eq_cbdc(1,1),max(p_grid))/y_func(eq_cbdc(1,2)+eq_cbdc(1,1)))/(alpha1+alpha2+alpha3);
    
Output_CBDC = alpha1*eq_cbdc(:,1)+alpha2*eliquidity+alpha3*min(eq_cbdc(:,1)+eliquidity,p_star)+B+eq_cbdc(:,end)+f(eq_cbdc(:,end))-eq_cbdc(:,2)...
    +(eq_cbdc(:,2)./(1+eq_cbdc(:,4))-eq_cbdc(:,end))*(1+i_reserve)/(1+pi_m);
Output_CBDCR = alpha1*eq_cbdc_R(:,1)+alpha2*eliquidityR+alpha3*min(eq_cbdc_R(:,1)+eliquidityR,p_star)+B+eq_cbdc_R(:,end)+f(eq_cbdc_R(:,end))-eq_cbdc_R(:,2)...
     +(eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4))-eq_cbdc_R(:,end)).*(1+max(i_reserve,icbdc_grid(:)))/(1+pi_m);
% CMoutput = eq_cbdc(:,end)+f(eq_cbdc(:,end))-eq_cbdc(:,2)+B;
% 
% CMoutputR = eq_cbdc_R(:,end)+f(eq_cbdc_R(:,end))-eq_cbdc_R(:,2)+B;
% id = find(CMoutput/CMoutput(1)>1,1,'last')
% icbdc_grid(id)
% figure
% hold on
% %plot(icbdc_grid*100,eliquidity)
% %plot(icbdc_grid*100, alpha1*eq_cbdc(:,1))
% %plot(icbdc_grid*100, alpha2*eliquidity)
% plot(icbdc_grid*100, alpha1*eq_cbdc(:,1)+alpha2*eliquidity+alpha3*min(eq_cbdc(:,1)+eliquidity,p_star))
% plot(icbdc_grid*100,eq_cbdc(:,end)+f(eq_cbdc(:,end))-eq_cbdc(:,2))





output_index=(Output_CBDC/Output_CBDC(1)-1)*100;
output_indexR=(Output_CBDCR/Output_CBDCR(1)-1)*100;

%% both as reserve and not as reserve
macro_implication=figure
psi_D = eq_cbdc(:,2)./(1+eq_cbdc(:,4));
psi_D_R = eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4));
subplot(3,3,1)
hold on
h1=plot(icbdc_grid*100,(psi_D./psi_D(1)-1)*100,'linewidth',2);
h2=plot(icbdc_grid*100,(psi_D_R./psi_D_R(1)-1)*100,'--red','linewidth',2);
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit')

subplot(3,3,2)
hold on
plot(icbdc_grid*100,(eq_cbdc(:,end)./eq_cbdc(1,end)-1)*100,'linewidth',2)
plot(icbdc_grid*100,(eq_cbdc_R(:,end)./eq_cbdc(1,end)-1)*100,'--red','linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan')

subplot(3,3,3)
hold on
plot(icbdc_grid*100,output_index,'linewidth',2)
plot(icbdc_grid*100,output_indexR,'--red','linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Output')

ndrate=((1+eq_cbdc(:,4))*(1+pi_m)-1)*100;
ndrateR=((1+eq_cbdc_R(:,4))*(1+pi_m)-1)*100;
subplot(3,3,4)
hold on
plot(icbdc_grid*100,((1+eq_cbdc(:,4))*(1+pi_m)-1)*100,'linewidth',2)
plot(icbdc_grid*100,((1+eq_cbdc_R(:,4))*(1+pi_m)-1)*100,'--red','linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit Rate')

nlrate=((1+eq_cbdc(:,3))*(1+pi_m)-1)*100;
nlrateR=((1+eq_cbdc_R(:,3))*(1+pi_m)-1)*100;
subplot(3,3,5)
hold on
plot(icbdc_grid*100,((1+eq_cbdc(:,3))*(1+pi_m)-1)*100,'linewidth',2)
plot(icbdc_grid*100,((1+eq_cbdc_R(:,3))*(1+pi_m)-1)*100,'--red','linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan Rate')

subplot(3,3,6)
hold on
plot(icbdc_grid*100,nlrate-ndrate,'linewidth',2)
plot(icbdc_grid*100,nlrateR-ndrateR,'--red','linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Spread')

hL = subplot(3,3,8);
poshL = get(hL,'position');     % Getting its position

lgd = legend(hL,[h1,h2],'CBDC not as Reserve','CBDC as Reserve', 'NumColumns',2);
set(lgd,'position',poshL);      % Adjusting legend's position
axis(hL,'off');                 % Turning its axis off

saveas(macro_implication,[location,'interest_bearing_',name],'epsc')
saveas(macro_implication,[location,'interest_bearing_',name],'png')

%% Not as reserves
macro_implication=figure
psi_D = eq_cbdc(:,2)./(1+eq_cbdc(:,4));
psi_D_R = eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4));
subplot(2,3,4)
hold on
h1=plot(icbdc_grid*100,(psi_D./psi_D(1)-1)*100,'linewidth',2);
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(2,3,5)
hold on
plot(icbdc_grid*100,(eq_cbdc(:,end)./eq_cbdc(1,end)-1)*100,'linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(2,3,6)
hold on
plot(icbdc_grid*100,output_index,'linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Output')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

ndrate=((1+eq_cbdc(:,4))*(1+pi_m)-1)*100;
ndrateR=((1+eq_cbdc_R(:,4))*(1+pi_m)-1)*100;
subplot(2,3,1)
hold on
plot(icbdc_grid*100,((1+eq_cbdc(:,4))*(1+pi_m)-1)*100,'linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Deposit Rate')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

nlrate=((1+eq_cbdc(:,3))*(1+pi_m)-1)*100;
nlrateR=((1+eq_cbdc_R(:,3))*(1+pi_m)-1)*100;
subplot(2,3,2)
hold on
plot(icbdc_grid*100,((1+eq_cbdc(:,3))*(1+pi_m)-1)*100,'linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Loan Rate')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

subplot(2,3,3)
hold on
plot(icbdc_grid*100,nlrate-ndrate,'linewidth',2)
xlim([-.5,2])
xlabel('$i_e$', 'interpreter','latex')
ylabel('Spread')
ax=gca
ax.XAxis.FontSize = 15;
ax.YAxis.FontSize = 15;

saveas(macro_implication,[location,'interest_bearing_nonreserve_',name],'epsc') % this line generates Figure 6 in the paper
saveas(macro_implication,[location,'interest_bearing_non_reserve_',name],'png')


welfare_func=@(x)alpha1*(u(x*y_func(eq_cbdc(1,1)))-p_func((eq_cbdc(1,1))))+alpha2*(u(x*y_func(eq_cbdc(1,2)))-p_func((eq_cbdc(1,2))))...
+alpha3*(u(x*y_func(eq_cbdc(1,1)+eq_cbdc(1,2)))-p_func((eq_cbdc(1,1)+eq_cbdc(1,2))))+B/2*log(x*B/2)-(B/2-eq_cbdc(1,2)+eq_cbdc(1,2)/(1+eq_cbdc(1,4)))...
-((1+i_reserve)/(1+pi_m)-1).*(eq_cbdc(1,2)./(1+eq_cbdc(1,4))-eq_cbdc(1,end));



%% buyers bear the operational costs of CBDC 


%alternative welfare calculation
welfare_cost=alpha1*(u(y_func(eq_cbdc(:,1)))-p_func((eq_cbdc(:,1))))+alpha2*(u(y_func(eliquidity))-p_func((eliquidity)))...
+alpha3*(u(y_func(eq_cbdc(:,1)+eliquidity))-p_func((eq_cbdc(:,1)+eliquidity)))-(B/2+eq_cbdc(:,2)./(1+eq_cbdc(:,4))/beta-eq_cbdc(:,2))+B/2*log(B/2)...
- (eliquidity - eq_cbdc(:,2))./(1+eq_cbdc(:,4))*cost;
welfareR_cost=alpha1*(u(y_func(eq_cbdc_R(:,1)))-p_func((eq_cbdc_R(:,1))))+alpha2*(u(y_func(eliquidityR))-p_func((eliquidityR)))...
+alpha3*(u(y_func(eq_cbdc_R(:,1)+eliquidityR))-p_func((eq_cbdc_R(:,1)+eliquidityR)))+B/2*log(B/2)-(B/2+eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4))/beta-eq_cbdc_R(:,2))...
- (eliquidityR - eq_cbdc_R(:,2))./(1+eq_cbdc_R(:,4))*cost;


for i=1:length(icbdc_grid)
    i
    welfare_comp(i) = fzero(@(x)welfare_func(x)-welfare_cost(i),[0.5,1.5]);
    welfare_compR(i) = fzero(@(x)welfare_func(x)-welfareR_cost(i),[0.5,1.5]);
end
h3=figure
hold on
plot(icbdc_grid*100,(welfare_comp-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_compR-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,5])
legend('CBCD not Reserve','CBDC as Reserve','location','northwest')
figure_halfpage
saveas(h3,[location,'welfare_b_cost',name],'png')




%entrepreneur welfare
welfare_e = f(eq_cbdc(:,end))-(1+eq_cbdc(:,3)).*eq_cbdc(:,end);

welfare_eR = f(eq_cbdc_R(:,end))-(1+eq_cbdc_R(:,3)).*eq_cbdc_R(:,end);

h2=figure
hold on
plot(icbdc_grid*100,(welfare_e/welfare_e(1)-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_eR/welfare_eR(1)-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,5])
[m,n]=max(welfare_e/welfare_e(1))
[m,nR]=max(welfare_eR/welfare_eR(1))
figure_halfpage
legend('CBDC not Reserve','CBDC as Reserve','location','southeast')

saveas(h2,[location,'welfare_e',name],'png')


%banker welfare
welfare_banker = max(1+eq_cbdc(:,3)-1/(1+pi_eff),0).*eq_cbdc(:,end)-eq_cbdc(:,2)-(cost-(1/(1+pi_eff))).*eq_cbdc(:,2)./(1+eq_cbdc(:,4));
welfare_bankerR = max(1+eq_cbdc_R(:,3)-max(1/(1+pi_eff),(1+icbdc_grid(:))/(1+pi_m)),0).*eq_cbdc_R(:,end)-eq_cbdc_R(:,2)-(cost-max(1/(1+pi_eff),(1+icbdc_grid(:))/(1+pi_m))).*eq_cbdc_R(:,2)./(1+eq_cbdc_R(:,4));
h1=figure
hold on
plot(icbdc_grid*100,(welfare_banker/welfare_banker(1)-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_bankerR/welfare_bankerR(1)-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,2])
figure_halfpage
legend('CBDC not Reserve','CBDC as Reserve')
saveas(h1,[location,'welfare_bankers',name],'png')

%welfare of sellers

welfareS=alpha1*(-c(y_func(eq_cbdc(:,1)))+p_func((eq_cbdc(:,1))))+alpha2*(-c(y_func(eliquidity))+p_func((eliquidity)))...
+alpha3*(-c(y_func(eq_cbdc(:,1)+eliquidity))+p_func((eq_cbdc(:,1)+eliquidity)))-B/2+B/2*log(B/2);
welfareSR=alpha1*(-c(y_func(eq_cbdc_R(:,1)))+p_func((eq_cbdc_R(:,1))))+alpha2*(-c(y_func(eliquidityR))+p_func((eliquidityR)))...
+alpha3*(-c(y_func(eq_cbdc_R(:,1)+eliquidityR))+p_func((eq_cbdc_R(:,1)+eliquidityR)))+B/2*log(B/2)-B/2;
welfare_funcS=@(x)alpha1*(-c(y_func(eq_cbdc(1,1)))+p_func((eq_cbdc(1,1))))+alpha2*(-c(y_func(eq_cbdc(1,2)))+p_func((eq_cbdc(1,2))))...
+alpha3*(-c(y_func(eq_cbdc(1,1)+eq_cbdc(1,2)))+p_func((eq_cbdc(1,1)+eq_cbdc(1,2))))+B/2*log(x*B/2)-B/2;
for i=1:length(icbdc_grid)
    i
    welfare_compS(i) = fzero(@(x)welfare_funcS(x)-welfareS(i),[0.5,1.5]);
    welfare_compSR(i) = fzero(@(x)welfare_funcS(x)-welfareSR(i),[0.5,1.5]);
end
h4=figure
hold on
plot(icbdc_grid*100,(welfare_compS-1)*100,'linewidth',2)
plot(icbdc_grid*100,(welfare_compSR-1)*100,'--red','linewidth',2)
xlabel('$i_e$','interpreter','latex')
xlim([0,5])
legend('CBCD not Reserve','CBDC as Reserve','location','northwest')
figure_halfpage
saveas(h4,[location,'welfare_S',name],'png')


%=============================
%Compute certain numbers
%=============================
[loan_max,nmax] = max(eq_cbdc(:,end)./eq_cbdc(1,end));
[loan_maxR,nmaxR] = max(eq_cbdc_R(:,end)./eq_cbdc(1,end));
i_max=icbdc_grid(nmax);
i_maxR=icbdc_grid(nmaxR);
id = find(eq_cbdc(:,end)./eq_cbdc(1,end)>=1,1,'last');
id1 = find(eq_cbdc(:,end)./eq_cbdc(1,end)>1+1e-8,1,'first');
idR1 = find(eq_cbdc_R(:,end)./eq_cbdc(1,end)>1,1,'first');

idR = find(eq_cbdc_R(:,end)./eq_cbdc(1,end)>=1,1,'last');
i_max1=icbdc_grid(id);
i_min1=icbdc_grid(id1);
i_minR1=icbdc_grid(idR1);
i_maxR1=icbdc_grid(idR);

[output_max,nmax] = max(output_index);
[outputR_max,nmaxR] = max(output_indexR);
ioutput_max=icbdc_grid(nmax);
ioutput_maxR=icbdc_grid(nmaxR);
id = find(output_index>=0,1,'last');
id1 = find(output_index>0,1,'first');
idR = find(output_indexR>=0,1,'last');
idR1 = find(output_indexR>0,1,'first');

iy_max1=icbdc_grid(id);
iy_min1=icbdc_grid(id1);
iy_maxR1=icbdc_grid(idR);
iy_minR1=icbdc_grid(idR1);


disp('=====================================================================================');
disp(['If CBDC does not serve as reserves, it increases Bank Lending if its rate is between 0.3049% and ' num2str(i_max1*100) '%']);
disp(['It increases output if its rate is between 0.3049% and ' num2str(iy_max1*100) '%']);
disp(['The CBDC rate that maximizes lending is ' num2str(i_max*100) '% and that maximizes output is ' num2str(ioutput_max*100) '%' ]);
disp(['The maximum increase in lending is ' num2str((loan_max-1)*100) '% and the maximum increase in output is ' num2str(output_max) '%' ]);
disp('')
disp('=====================================================================================');
disp(['If CBDC serves as reserves, it increases Bank Lending if its rate is between 0.3049% and ' num2str(i_maxR1*100) '%']);
disp(['It increases output if its rate is between 0.3049% and ' num2str(iy_maxR1*100) '%']);
disp(['The CBDC rate that maximizes lending is ' num2str(i_maxR*100) '% and that maximizes output is ' num2str(ioutput_maxR*100) '%' ]);
disp(['The maximum increase in lending is ' num2str((loan_maxR-1)*100) '% and the maximum increase in output is ' num2str(outputR_max) '%' ]);

